Dissimilarity time series analysis

Primary script 4a for Kitchel et al. 2023 in prep taxonomic diversity manuscript.

library(data.table)
library(vegan)
library(sf)
library(concaveman) #polygon around points
library(betapart) #allows us to partition beta diversity
library(geosphere)
library(ggpubr) #stat_regline_equation
library(nlme)
library(mgcv) #to make gam
library(cowplot)
library(lme4)
library(stringr)

Pull dissimilarity table

#Pull Dissimilarity Means
distances_dissimilarities_allyears.r <- readRDS(here::here("output", "dissimilarities", "distances_dissimilarities_allyears.r.rds"))

#make survey and survey unit factors
distances_dissimilarities_allyears.r[,survey:=factor(survey)][,survey_unit:=factor(survey_unit)]

#adjust years
distances_dissimilarities_allyears.r[,year_adj := year-min(year)+1]

#add new variable for year in sequence per region
distances_dissimilarities_allyears.r[,first_year := min(year),.(survey_unit)]
distances_dissimilarities_allyears.r[,last_year := max(year),.(survey_unit)]

#distances_dissimilarities_allyears.r[,year_in_seq := year-first_year+1]

distances_dissimilarities_allyears.r[,years_sampled := last_year-first_year+1]

###Palette for Plotting Palette for plotting all 42 survey units

survey_unit.list <- levels(factor(distances_dissimilarities_allyears.r[,survey_unit]))

palette_42 <- c(
  "#5A5156", #AI
  "#DF00DB", #BITS-1
  "#DB8EDA", #BITS-4
  "#F6222E", #CHL
  "#F8A19F", #DFO-NF
  "#16FF32", #DFO-QCS
  "#325A9B", #EBS
  "#3283FE", #EVHOE
  "#FEAF16", #FR-CGFS
  "#fccb6d", #GMEX-Fall
  "#1C8356", #GMEX-Summer
  "#C4451C", #GOA
  "#85660D", #GRL-DE
  "#B0009F", #GSL-N
  "#BF79B8", #GSL-S
  "#1CBE4F", #ICE-GFS
  "#782AB6", #IE-IGFS
  "#90AD1C", #MEDITS
  "#6B003A", #NAM
  "#A75B00", #NEUS-Fall
  "#E3B072", #NEUS-Spring
  "#02E8B6", #NIGFS-1
  "#97E7D5", #NIGFS-4
  "#B00068", #Nor-BTS-3
  "#00B9E3", #NS-IBTS-1
  "#95E2F4", #NS-IBTS-3
  "#B3CE73", #NZ-CHAT
  "#689500", #NZ-ECSI
  "#364d02",#NZ-SUBA
  "#AAF400", #NZ-WCSI
  "#AA0DFE", #PT-IBTS
  "#7f9eb8", #ROCKALL
  "#FA0087", #S-GEORG
  "#DEA0FD", #SCS-Summer
  "#FCEF88", #SEUS-fall
  "#A59405", #SEUS-spring
  "#FCE100", #SEUS-summer
  "#544563", #SWC-IBTS-1
  "#a37fc7", #SWC-IBTS-4
  "#C075A6", #WCANN
  "#BDCDFF", #ZAF-ATL
  "#003EFF"  #ZAF-IND
)

color_link <- data.table(survey_unit = survey_unit.list,hex = palette_42)

Add names for plotting

name_helper <- data.table(Survey_Name_Season = 
                            c(
                              "Aleutian Islands",
                                    "Baltic Sea Q1",
                                    "Baltic Sea Q4",
                                    "Chile",
                                    "Newfoundland",
                                    "Queen Charlotte Sound",
                                    "Eastern Bering Sea",
                                    "Bay of Biscay",
                                    "English Channel",
                                    "Gulf of Mexico Summer",
                                    "Gulf of Alaska",
                                    "Greenland",
                                    "N Gulf of St. Lawrence",
                                    "S Gulf of St. Lawrence",
                                    "Iceland",
                                    "Irish Sea",
                                    "Mediterranean",
                                    "Namibia",
                                    "NE US Fall",
                                    "NE US Spring",
                                    "N Ireland Q1",
                                    "N Ireland Q4",
                                    "Barents Sea Norway Q3",
                                    "N Sea Q1",
                                    "N Sea Q3",
                                    "Chatham Rise NZ",
                                    "E Coast S Island NZ",
                                    "W Coast S Island NZ",
                                    "Portugal",
                                    "S Georgia",
                                  "Scotian Shelf Summer",
                                  "SE US Fall",
                                  "SE US Spring",
                                  "SE US Summer",
                                  "W Coast US",
                                  "Atlantic Ocean ZA",
                                  "Indian Ocean ZA",
                                   "Rockall Plateau",
                                  "Scotland Shelf Sea Q1",
                                  "Scotland Shelf Sea Q4",
                                  "Falkland Islands",
                                  "Gulf of Mexico Fall",
                                  "Barents Sea Norway Q1",
                                  "Sub-Antarctic NZ",
                                  "Scotian Shelf Spring"),
                          survey_unit = c(
                                  "AI",        
                                  "BITS-1",    
                                  "BITS-4",    
                                  "CHL",       
                                  "DFO-NF",    
                                  "DFO-QCS",   
                                  "EBS",       
                                  "EVHOE",     
                                  "FR-CGFS",   
                                  "GMEX-Summer",
                                  "GOA",       
                                  "GRL-DE",    
                                  "GSL-N",     
                                  "GSL-S",     
                                  "ICE-GFS",   
                                  "IE-IGFS",   
                                  "MEDITS",    
                                  "NAM",       
                                  "NEUS-Fall", 
                                  "NEUS-Spring",
                                  "NIGFS-1",   
                                  "NIGFS-4",   
                                  "Nor-BTS-3", 
                                  "NS-IBTS-1", 
                                  "NS-IBTS-3", 
                                  "NZ-CHAT",   
                                  "NZ-ECSI",   
                                  "NZ-WCSI",   
                                  "PT-IBTS",   
                                  "S-GEORG",   
                                  "SCS-SUMMER",
                                  "SEUS-fall", 
                                  "SEUS-spring",
                                  "SEUS-summer",
                                  "WCANN",     
                                  "ZAF-ATL",   
                                  "ZAF-IND",
                                  "ROCKALL",
                                  "SWC-IBTS-1",
                                  "SWC-IBTS-4",
                                  "FALK",
                                  "GMEX-Fall",
                                  "Nor-BTS-1",
                                  "NZ-SUBA",
                                  "SCS-SPRING"
                          ))


color_link <- name_helper[color_link, on = "survey_unit"]

#save to pull into other scripts
fwrite(color_link, here::here("output","color_link.csv"))

##Make GAMs

Bray Curtis

bray_curtis_total_gam <- gam(bray_curtis_dissimilarity_total_mean ~ year + s(survey_unit, year, bs = "fs", m = 1),#random smooth
                            data = distances_dissimilarities_allyears.r)

Jaccard

jaccard_total_gam <- gam(jaccard_dissimilarity_total_mean ~ year + s(survey_unit, year, bs = "fs", m = 1),#random smooth
                            data = distances_dissimilarities_allyears.r)

##Make LMERS

Bray These all converged

#running with lme instead of lmer gave same results, but allowed for calculation of p-value
bray_curtis_total_lme <- lme(bray_curtis_dissimilarity_total_mean ~ year_adj, random = (~1 + year_adj|survey_unit),data = distances_dissimilarities_allyears.r)

#but also run with lmer for confint
#total
bray_curtis_total_lmer <- lmer(bray_curtis_dissimilarity_total_mean ~ year_adj + (1 + year_adj|survey_unit),data = distances_dissimilarities_allyears.r)

#balanced (for comparison plot)
bray_curtis_balanced_lmer <- lmer(bray_curtis_dissimilarity_balanced_mean ~ year_adj + (1 + year_adj|survey_unit),data = distances_dissimilarities_allyears.r)

Warning: Model failed to converge with max|grad| = 0.0930452 (tol = 0.002, component 1)

#gradient (for comparison plot)
bray_curtis_gradient_lmer <- lmer(bray_curtis_dissimilarity_gradient_mean ~ year_adj + (1 + year_adj|survey_unit),data = distances_dissimilarities_allyears.r)

Warning: Model failed to converge with max|grad| = 0.5086 (tol = 0.002, component 1)

#Jaccard total (for comparison plot)
jaccard_total_lmer <- lmer(jaccard_dissimilarity_total_mean ~ year_adj + (1 + year_adj|survey_unit),data = distances_dissimilarities_allyears.r)

Warning: Model failed to converge with max|grad| = 0.163399 (tol = 0.002, component 1)

summary(bray_curtis_total_lme)

Linear mixed-effects model fit by REML
  Data: distances_dissimilarities_allyears.r 

Random effects:
 Formula: ~1 + year_adj | survey_unit
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev      Corr  
(Intercept) 0.106024695 (Intr)
year_adj    0.001955226 -0.769
Residual    0.024352895       

Fixed effects:  bray_curtis_dissimilarity_total_mean ~ year_adj 

 Correlation: 
         (Intr)
year_adj -0.798

Standardized Within-Group Residuals:
        Min          Q1         Med          Q3         Max 
-4.93606507 -0.46317561  0.02981855  0.50587465  4.55375745 

Number of Observations: 896
Number of Groups: 42 

anova(bray_curtis_total_lme)

bray_curtis_total_coefs <- data.table(coef(bray_curtis_total_lme))
bray_curtis_total_coefs[,survey_unit := rownames(coef(bray_curtis_total_lme))][,Year := round(year_adj,5)][,Intercept := round(`(Intercept)`,2)]
#View(bray_curtis_total_coefs)

bray_curtis_total_coefs <- bray_curtis_total_coefs[color_link, on = "survey_unit"]

bray_curtis_total_coefs.exp <- bray_curtis_total_coefs[,.(Survey_Name_Season, Intercept, Year)]

#export this table
fwrite(bray_curtis_total_coefs.exp, file = here::here("output","bray_curtis_total_coefs.exp.csv"))

Also, build a simple linear model with an interaction (instead of random slope and intercept for survey) THIS IS HOW WE CLASSIFY SURVEYS

bray_curtis_total_lm <- lm(bray_curtis_dissimilarity_total_mean ~ year_adj*survey_unit,data = distances_dissimilarities_allyears.r)

bray_curtis_total_coefs <- data.table(summary(bray_curtis_total_lm)$coefficients)
  bray_curtis_total_coefs[,var := rownames(summary(bray_curtis_total_lm)$coefficients)]
  
  #limit to interactions only (check this if there are any model changes!) row 2 and rows 44:84
  bray_curtis_total_coefs.r <- bray_curtis_total_coefs[c(2,44:84),]
  
  #adjust survey unit name by deleting beginning of string
  bray_curtis_total_coefs.r[,survey_unit := substr(var, 21, str_length(var))][var == "year_adj",survey_unit := "AI"]
  
  #calculate interaction coefficients
  AI_estimate <- bray_curtis_total_coefs.r[1,Estimate]
  bray_curtis_total_coefs.r[1,survey_unit_coefficient := AI_estimate]
  bray_curtis_total_coefs.r[2:42,survey_unit_coefficient := (AI_estimate + Estimate)]
  
  #homogenizing vs differentiating
  bray_curtis_total_coefs.r[,differentiating := ifelse((`Std. Error`< abs(survey_unit_coefficient) & survey_unit_coefficient > 0),1,0)][,homogenizing := ifelse((`Std. Error`< abs(survey_unit_coefficient) & survey_unit_coefficient < 0),1,0)]
  
  #lower and upper CI
    bray_curtis_total_coefs.r[,lwr := survey_unit_coefficient-`Std. Error`][,upr:= survey_unit_coefficient+`Std. Error`]
  
  bray_curtis_total_coefs_sum <- data.table(differentiating = sum(  bray_curtis_total_coefs.r$differentiating), homogenizing = sum(  bray_curtis_total_coefs.r$homogenizing))
  
  #hex by trend
  bray_curtis_total_coefs.r[,trend_color := ifelse(homogenizing == 1, "#e7ac5b", ifelse(differentiating == 1,"#91c874","#cbbfde"))]
  
  bray_curtis_total_coefs.r_alpha <- setorder(bray_curtis_total_coefs.r, Survey_Name_Season)

Error in setorderv(x, cols, order, na.last) : 
  some columns are not in the data.table: Survey_Name_Season

Get LMER model as predictions

# need to sort out year in seq versus overall year models
#new data for lmer
lmer_year <- seq(min(distances_dissimilarities_allyears.r[,year]), max(distances_dissimilarities_allyears.r[,year]), by = 0.1)

lmer_year_adj <- seq(min(distances_dissimilarities_allyears.r[,year_adj]), max(distances_dissimilarities_allyears.r[,year_adj]), by = 0.1)

#predict average lmer
lmer_bray_total_predictions <- data.table(year = lmer_year, year_adj = lmer_year_adj)

#confidence intervals
bray_curtis_total_lmer_confint <- confint(bray_curtis_total_lmer)

Computing profile confidence intervals ...

#populate data table of lmer predictions
lmer_bray_total_predictions[,bray_curtis_lmer_preds := fixef(bray_curtis_total_lmer)[[1]] + fixef(bray_curtis_total_lmer)[[2]] * year_adj][,bray_curtis_lmer_preds_lowerCI := bray_curtis_total_lmer_confint[5] + bray_curtis_total_lmer_confint[6] * year_adj][,bray_curtis_lmer_preds_upperCI := bray_curtis_total_lmer_confint[11] + bray_curtis_total_lmer_confint[12] * year_adj]

# need to sort out year in seq versus overall year models
#new data for lmer


#predict average lmer
lmer_bray_balanced_predictions <- data.table(year = lmer_year, year_adj = lmer_year_adj)

#confidence intervals
bray_curtis_balanced_lmer_confint <- confint(bray_curtis_balanced_lmer)

Computing profile confidence intervals ...

#populate data table of lmer predictions
lmer_bray_balanced_predictions[,bray_curtis_lmer_preds := fixef(bray_curtis_balanced_lmer)[[1]] + fixef(bray_curtis_balanced_lmer)[[2]] * year_adj][,bray_curtis_lmer_preds_lowerCI := bray_curtis_balanced_lmer_confint[5] + bray_curtis_balanced_lmer_confint[6] * year_adj][,bray_curtis_lmer_preds_upperCI := bray_curtis_balanced_lmer_confint[11] + bray_curtis_balanced_lmer_confint[12] * year_adj]

# need to sort out year in seq versus overall year models


#predict average lmer
lmer_bray_gradient_predictions <- data.table(year = lmer_year, year_adj = lmer_year_adj)

#confidence intervals
bray_curtis_gradient_lmer_confint <- confint(bray_curtis_gradient_lmer)

Computing profile confidence intervals ...

#populate data table of lmer predictions
lmer_bray_gradient_predictions[,bray_curtis_lmer_preds := fixef(bray_curtis_gradient_lmer)[[1]] + fixef(bray_curtis_gradient_lmer)[[2]] * year_adj][,bray_curtis_lmer_preds_lowerCI := bray_curtis_gradient_lmer_confint[5] + bray_curtis_gradient_lmer_confint[6] * year_adj][,bray_curtis_lmer_preds_upperCI := bray_curtis_gradient_lmer_confint[11] + bray_curtis_gradient_lmer_confint[12] * year_adj]

# need to sort out year in seq versus overall year models


#predict average lmer
lmer_jaccard_total_predictions <- data.table(year = lmer_year, year_adj = lmer_year_adj)

#confidence intervals
jaccard_total_lmer_confint <- confint(jaccard_total_lmer)

Computing profile confidence intervals ...

#populate data table of lmer predictions
lmer_jaccard_total_predictions[,jaccard_lmer_preds := fixef(jaccard_total_lmer)[[1]] + fixef(jaccard_total_lmer)[[2]] * year_adj][,jaccard_lmer_preds_lowerCI := jaccard_total_lmer_confint[5] + jaccard_total_lmer_confint[6] * year_adj][,jaccard_lmer_preds_upperCI := jaccard_total_lmer_confint[11] + jaccard_total_lmer_confint[12] * year_adj]

###Move forward with Bray Curtis total (others to supplement, although will include in plot below) Coefficients for LMER by survey_unit

#unique survey unit year combos
survey_unit_sampling_years <- unique(distances_dissimilarities_allyears.r[,.(survey_unit, year_adj, year, years_sampled)])

#just # years
survey_unit_years <- unique(survey_unit_sampling_years[,.(survey_unit, years_sampled)])

# see group coefficients
model_coefs_reduced <- data.table(transform(as.data.frame(ranef(bray_curtis_total_lmer)), lwr = condval - 1.96*condsd, upr = condval + 1.96*condsd))
#https://stackoverflow.com/questions/69805532/extract-the-confidence-intervals-of-lmer-random-effects-plotted-with-dotplotra


#ONLY SLOPES
model_coefs_reduced <- model_coefs_reduced[term == "year_adj",]

model_coefs_reduced[,survey_unit := grp][,year_adj := condval]

#merge with duration of survey
model_coefs_reduced_length <- model_coefs_reduced[survey_unit_sampling_years, on = "survey_unit"]


model_coefs_reduced_length[,years_sampled := as.numeric(years_sampled)][,Directional_Change := ifelse(year_adj > 0, "Differentiation","Homogenization")]

#does it cross zero?
model_coefs_reduced_length[,significant := ifelse(lwr >0 & upr>0,T,ifelse(lwr<0 & upr<0,T,F))]

#significant directional change
model_coefs_reduced_length[,Directional_Change_sig := ifelse(year_adj > 0 & significant == T, "Differentiation",ifelse(year_adj < 0 & significant == T, "Homogenization", "No trend in\ndissimilarity"))]


#min max distances_dissimilarities
min_bray_reduced <- min(distances_dissimilarities_allyears.r[,bray_curtis_dissimilarity_total_mean], na.rm = T)
max_bray_reduced <- max(distances_dissimilarities_allyears.r[,bray_curtis_dissimilarity_total_mean], na.rm = T)

model_coefs_reduced_length <- model_coefs_reduced_length[color_link, on = "survey_unit"]

#delete any NAs
model_coefs_reduced_length <- na.omit(model_coefs_reduced_length, cols = "significant")

#order table by coefficient
setorder(model_coefs_reduced_length, year_adj)

BC_total_model_coefs_reduced_length.unique <- unique(model_coefs_reduced_length[,.(condval,condsd, lwr, upr, survey_unit, year_adj, years_sampled, Directional_Change, hex, Survey_Name_Season, significant, Directional_Change_sig)]) 

#extract color hexes
#year adj coef order
color_year_adj_order <- BC_total_model_coefs_reduced_length.unique[,hex]

#alphabetical order
BC_total_model_coefs_reduced_length.unique.alpha <- setorder(BC_total_model_coefs_reduced_length.unique, Survey_Name_Season)

BC_total_model_coefs_reduced_length.unique.alpha[,trend_color := ifelse(Directional_Change_sig == "Homogenization", "#e7ac5b", ifelse(Directional_Change_sig == "Differentiation","#91c874","#cbbfde"))]

color_alpha_order <- BC_total_model_coefs_reduced_length.unique.alpha[,hex]
color_alpha_order_bytrend <- BC_total_model_coefs_reduced_length.unique.alpha[, trend_color]

saveRDS(BC_total_model_coefs_reduced_length.unique, here::here("output","region_stats","BC_total_model_coefs_reduced_length.unique.Rds"))

Coefficients by interactions of linear model

#extract color hexes
#year adj coef order
color_year_adj_order <- BC_total_model_coefs_reduced_length.unique[,hex]

#alphabetical order
BC_total_model_coefs_reduced_length.unique.alpha <- setorder(BC_total_model_coefs_reduced_length.unique, Survey_Name_Season)

BC_total_model_coefs_reduced_length.unique.alpha[,trend_color := ifelse(Directional_Change_sig == "Homogenization", "#e7ac5b", ifelse(Directional_Change_sig == "Differentiation","#91c874","#cbbfde"))]

color_alpha_order <- BC_total_model_coefs_reduced_length.unique.alpha[,hex]
color_alpha_order_bytrend <- BC_total_model_coefs_reduced_length.unique.alpha[, trend_color]

#key trend color and survey
trend_color_survey <- unique(BC_total_model_coefs_reduced_length.unique.alpha[,.(survey_unit, Survey_Name_Season,hex, trend_color)])

saveRDS(BC_total_model_coefs_reduced_length.unique, here::here("output","region_stats","BC_total_model_coefs_reduced_length.unique.Rds"))

Total versus balanced BC plot

ggplot(distances_dissimilarities_allyears.r) +
  geom_point(aes(bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_balanced_mean)) +
  geom_abline(slope = 1, intercpet = 0) +
  geom_smooth(aes(bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_balanced_mean)) +
  theme_classic() +
  labs(x = "Total BC dissimilarity",  y = "Balanced changes in abundance/biomass")

Warning: Ignoring unknown parameters: `intercpet`

ggplot(distances_dissimilarities_allyears.r) +
  geom_point(aes(bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_gradient_mean)) +
  geom_abline(slope = 1, intercpet = 0) +
  geom_smooth(aes(bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_gradient_mean)) +
  theme_classic() +
  labs(x = "Total BC dissimilarity",  y = "Abundance/biomass gradient")

Warning: Ignoring unknown parameters: `intercpet`

#all on one plot
#reduced
distances_dissimilarities_allyears.r.r <- distances_dissimilarities_allyears.r[,.(year, bray_curtis_dissimilarity_total_mean, bray_curtis_dissimilarity_gradient_mean, bray_curtis_dissimilarity_balanced_mean, jaccard_dissimilarity_total_mean)]

#wide to long
dissim_long <- melt(distances_dissimilarities_allyears.r.r, id.vars = c("year"), variable.name = "Dissimilarity metric")

dissim_long[,`Dissimilarity metric` := factor(`Dissimilarity metric`, levels = c( "bray_curtis_dissimilarity_total_mean","bray_curtis_dissimilarity_balanced_mean","bray_curtis_dissimilarity_gradient_mean", "jaccard_dissimilarity_total_mean"))]

BC_dissimilarity_components <- ggplot() +
  geom_point(data = dissim_long[`Dissimilarity metric` != "jaccard_dissimilarity_total_mean",], aes(x = year, y = value, color = `Dissimilarity metric`), alpha = 0.2, size = 1) +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "black", alpha = 0.3) +
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
  #balanced
    geom_ribbon(data = lmer_bray_gradient_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "cadetblue3", alpha = 0.3) +
  geom_line(data = lmer_bray_gradient_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "cadetblue3") +
  #gradient
    geom_ribbon(data = lmer_bray_balanced_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "darksalmon", alpha = 0.3) +
  geom_line(data = lmer_bray_balanced_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "darksalmon") +
  #fix. colors of points
  scale_color_manual(values = c("black","darksalmon","cadetblue3"), labels = c("Total","Balanced variation\nin abundance","Biomass gradient")) +
  labs(x = "Year", y = "Value") +
  lims(y = c(min(dissim_long$value), max(dissim_long$value))) +
  theme_classic() +
  theme(legend.position = c(0.17,0.37))

ggsave(BC_dissimilarity_components, path = here::here("figures"), filename=("BC_dissimilarity_components.jpg"), width = 6.5, height = 5, units = "in")

#Jaccard
Jaccard_dissimilarity <- ggplot() +
  geom_point(data = dissim_long[`Dissimilarity metric` == "jaccard_dissimilarity_total_mean"], aes(x = year, y = value), alpha = 0.2, size = 1) +
  geom_ribbon(data = lmer_jaccard_total_predictions, aes(x = year, ymin = jaccard_lmer_preds_lowerCI, ymax = jaccard_lmer_preds_upperCI), fill = "black", alpha = 0.3) +
  geom_line(data = lmer_jaccard_total_predictions, aes(x = year, y = jaccard_lmer_preds), color = "black") +
  labs(x = "Year", y = "Total Jaccard dissimilarity") +
  lims(y = c(min(dissim_long$value), max(dissim_long$value))) +
  theme_classic()

ggsave(Jaccard_dissimilarity, path = here::here("figures"), filename=("Jaccard_dissimilarity.jpg"), width = 6.5, height = 5, units = "in")

#merge
library(cowplot)
all_dissimilarities_over_time <- plot_grid(BC_dissimilarity_components + theme(legend.background = element_rect(fill = "transparent")),
                                           Jaccard_dissimilarity,
                                           ncol = 2, align = "hv", labels = c("a.","b."))

ggsave(all_dissimilarities_over_time, filename = "all_dissimilarities_over_time.jpg", path = here::here("figures"), width = 10, height = 4.5, units = "in")

NA
NA

Lollipop SE plot by linear trend experienced

#"#73BA4D","#E0962C","#cbbfde"

BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend <- ggplot() +
    geom_errorbar(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, ymin = lwr, ymax = upr), fill = "grey", width = 0) + #add confidence intervals
  geom_point(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, size = years_sampled, color = Directional_Change_sig), stat = 'identity') +
  scale_color_manual(values = c("#73BA4D","#E0962C","#cbbfde"), name = "Dissimilarity trend", guide="none") +
  scale_size_binned(range = c(1,8), name = "Survey period length\n(years)") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = seq(-0.005, 0.0075, by = 0.0025), labels = c("-0.005","","0", "", "0.005",  "")) +
  xlab("Survey unit") +
  ylab("β-diversity trend") + #total BC dissimilarity trend
  coord_flip() +
  theme_classic() +
  theme(axis.text.y = element_text(face = "bold"), axis.title.y = element_blank(), axis.text.x = element_text(size = 15), axis.title.x = element_text(size = 15), legend.position = c(0.3,0.8), legend.direction = "vertical", legend.text = element_text(size = 15), legend.title = element_text(size = 16))

Warning: Ignoring unknown parameters: `fill`Warning: Ignoring unknown aesthetics: labelWarning: Ignoring unknown aesthetics: label

directional_change_legend_plot_colorbytrend <- BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend + 
  theme(legend.position = "right", legend.background = element_rect(fill= "transparent"), 
         legend.text = element_text(size = 15), legend.title = element_text(size = 16)) +
  guides(colour = guide_legend(override.aes = list(size=6)), size = "none")

###############
#FOR LINEAR MODEL INSTEAD OF LINEAR MIXED EFFECTS
###############
bray_curtis_total_coefs.r <- bray_curtis_total_coefs.r[color_link, on = "survey_unit"]
bray_curtis_total_coefs.r <- survey_unit_years[bray_curtis_total_coefs.r, on="survey_unit"]

bray_curtis_total_coefs.r[,Directional_Change_sig := ifelse(differentiating > 0, "Differentiation",ifelse(homogenizing > 0 , "Homogenization", "No trend in\ndissimilarity"))]

(BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend_LINEAR_MODEL <- ggplot() +
    geom_errorbar(data = bray_curtis_total_coefs.r, aes(x = reorder(Survey_Name_Season, survey_unit_coefficient) , y = survey_unit_coefficient, label = Survey_Name_Season, ymin = lwr, ymax = upr), fill = "grey", width = 0) + #add confidence intervals
  geom_point(data = bray_curtis_total_coefs.r, aes(x = reorder(Survey_Name_Season, survey_unit_coefficient) , y = survey_unit_coefficient, label = Survey_Name_Season, size = years_sampled, color = Directional_Change_sig), stat = 'identity') +
  scale_color_manual(values = c("#73BA4D","#E0962C","#cbbfde"), name = "Dissimilarity trend", guide="none") +
  scale_size_binned(range = c(1,8), name = "Survey period length\n(years)") +
  geom_hline(yintercept = 0) +
  #scale_y_continuous(breaks = seq(-0.005, 0.0075, by = 0.0025), labels = c("-0.005","","0", "", "0.005",  "")) +
  xlab("Survey unit") +
  ylab("β-diversity trend") + #total BC dissimilarity trend
  coord_flip() +
  theme_classic() +
  theme(axis.text.y = element_text(face = "bold"), axis.title.y = element_blank(), axis.text.x = element_text(size = 15), axis.title.x = element_text(size = 15), legend.position = c(0.3,0.8), legend.direction = "vertical", legend.text = element_text(size = 15), legend.title = element_text(size = 16)))

Warning: Ignoring unknown parameters: `fill`Warning: Ignoring unknown aesthetics: labelWarning: Ignoring unknown aesthetics: label

ggsave(BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend_LINEAR_MODEL, 
       path = here::here("figures"), filename = "BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend_LINEAR_MODEL.jpg", height = 7, width = 7)

#Still need to edit
directional_change_legend_plot_colorbytrend <- BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend + 
  theme(legend.position = "right", legend.background = element_rect(fill= "transparent"), 
         legend.text = element_text(size = 15), legend.title = element_text(size = 16)) +
  guides(colour = guide_legend(override.aes = list(size=6)), size = "none")

Wavy Line Plot for GAMs

Generate predicted values (this datatable will hold both jaccard and bray curtis)

#add colors and names to full dissimilarity data table
distances_dissimilarities_allyears.r <- distances_dissimilarities_allyears.r[color_link, on = "survey_unit"]

#generate new data to smooth lines (need year and season survey combinations)
year_survey_unit_expand.dt <- data.table(survey_unit = as.character(NULL), year = as.numeric(NULL), year_adj = as.numeric(NULL ))

for (i in 1:length(survey_unit.list)) {
  #generate year vectors
  survey_unit_years <- unique(distances_dissimilarities_allyears.r[survey_unit == survey_unit.list[i],.(survey_unit, year, year_adj)])
  
  years <- seq(min(survey_unit_years[,year]), max(survey_unit_years[,year]), by = 0.1)
  
  year_adjust <- seq(min(survey_unit_years[,year_adj]), max(survey_unit_years[,year_adj]), by = 0.1)
  
  year_survey_unit_expand.dt_addition <- data.table(survey_unit = survey_unit.list[i], year = years, year_adj = year_adjust)
  
  year_survey_unit_expand.dt <- rbind(year_survey_unit_expand.dt, year_survey_unit_expand.dt_addition)
}

#add colors and names to full year and survey unit combination table
year_survey_unit_expand.dt <- year_survey_unit_expand.dt[color_link, on = "survey_unit"]

Get model as predictions Bray Curtis

#for plotting, get model as predictions
bray_curtis_total_gam_predictions <- predict(bray_curtis_total_gam, se.fit = TRUE, newdata = year_survey_unit_expand.dt)

#merge into table
year_survey_unit_expand.dt[,bray_glm_mod_fit := bray_curtis_total_gam_predictions$fit][,bray_glm_mod_fit_SE := bray_curtis_total_gam_predictions$se.fit]

Jaccard Get model as predictions

#for plotting, get model as predictions
jaccard_total_gam_predictions <- predict(jaccard_total_gam, se.fit = TRUE, newdata = year_survey_unit_expand.dt)

#merge into table
year_survey_unit_expand.dt[,jaccard_glm_mod_fit := jaccard_total_gam_predictions$fit][,jaccard_glm_mod_fit_SE := jaccard_total_gam_predictions$se.fit]

Produce Plot of GAM and mean LMER line

points_wavylines_bray_total_year_reduced_gam_nolmer <- ggplot() +
 # geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.2) +
  geom_point(data = na.omit(distances_dissimilarities_allyears.r,cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 color = Survey_Name_Season), alpha = 0.5, size = 1) +
    geom_line(data = na.omit(year_survey_unit_expand.dt,cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season)) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt,cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
 # geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =  color_alpha_order, name = "Survey Unit") +
  scale_fill_manual(values =  color_alpha_order, guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year])),
       y = c(0.15,0.9)) +
  xlab("Year") +
ylab("total BC dissimilarity") +
  theme(legend.position = "null")

points_wavylines_bray_total_year_reduced_gam_nolmer

ggsave(points_wavylines_bray_total_year_reduced_gam_nolmer, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_nolmer.jpg", height = 5, width = 5, unit = "in")

#with lmer

points_wavylines_bray_total_year_reduced_gam <- ggplot() +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.3) +
  geom_point(data = na.omit(distances_dissimilarities_allyears.r, cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =  color_alpha_order, name = "Survey Unit") +
  scale_fill_manual(values =  color_alpha_order, guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year]))) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam

ggsave(points_wavylines_bray_total_year_reduced_gam, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam.jpg", height = 6, width = 6, unit = "in")

#ALT
#plot all, but same color scheme (grey)
points_wavylines_bray_total_year_reduced_gam_greyscale <- ggplot() +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.3) +
  geom_point(data = distances_dissimilarities_allyears.r,
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = year_survey_unit_expand.dt,
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = year_survey_unit_expand.dt, aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =  rep("black", times = length(unique(distances_dissimilarities_allyears.r$Survey_Name_Season))), name = "Survey Unit") +
  scale_fill_manual(values =  rep("black", times = length(unique(distances_dissimilarities_allyears.r$Survey_Name_Season))), guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year]))
       ) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam_greyscale

ggsave(points_wavylines_bray_total_year_reduced_gam_greyscale, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_greyscale.jpg", height = 6, width = 6, unit = "in")

Alternative, color by trend

points_wavylines_bray_total_year_reduced_gam_colorbytrend <- ggplot() +
  geom_ribbon(data = lmer_bray_total_predictions, aes(x = year, ymin = bray_curtis_lmer_preds_lowerCI, ymax = bray_curtis_lmer_preds_upperCI), fill = "grey", alpha = 0.3) +
  geom_point(data = na.omit(distances_dissimilarities_allyears.r, cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
  geom_line(data = lmer_bray_total_predictions, aes(x = year, y = bray_curtis_lmer_preds), color = "black") +
    scale_color_manual(values =    linear_model_color_order_factorlevels$trend_color, name = "Survey Unit") +
  scale_fill_manual(values =    linear_model_color_order_factorlevels$trend_color, guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year])),y = c(0.49,0.95)) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam_colorbytrend

ggsave(points_wavylines_bray_total_year_reduced_gam_colorbytrend, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_colorbytrend.jpg", height = 6, width = 11, unit = "in")

points_wavylines_bray_total_year_reduced_gam_colorbytrend_squigglesonly <- ggplot() +
  geom_point(data = na.omit(distances_dissimilarities_allyears.r, cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), alpha = 0.4, size = 1, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season), alpha = 0.6) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt, cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
    scale_color_manual(values =    linear_model_color_order_factorlevels$trend_color, name = "Survey Unit") +
  scale_fill_manual(values =    linear_model_color_order_factorlevels$trend_color, guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year])), y = c(0.49,0.95)) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam_colorbytrend_squigglesonly

ggsave(points_wavylines_bray_total_year_reduced_gam_colorbytrend_squigglesonly, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_colorbytrend_squigglesonly.jpg", height = 6, width = 11, unit = "in")

BC

#plot each independently for supplement
#all survey names = 
all_survey_names <- sort(unique(color_link$Survey_Name_Season))
#list of plots
points_wavylines_bray_total_year_reduced_gam_individual <- list()
for (i in 1:length(all_survey_names)) {
points_wavylines_bray_total_year_reduced_gam_individual[[i]] <- ggplot() +
  geom_point(data = distances_dissimilarities_allyears.r[Survey_Name_Season == all_survey_names[i]],
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean), alpha = 0.4, color = "black") +
    geom_line(data = year_survey_unit_expand.dt[Survey_Name_Season == all_survey_names[i]],
             aes(x = year,
                 y = bray_glm_mod_fit), alpha = 0.6) +
  geom_ribbon(data = year_survey_unit_expand.dt[Survey_Name_Season == all_survey_names[i]], aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE), alpha=0.1) + #add standard error
  theme_classic() +
#  lims(x = c(min(distances_dissimilarities_allyears.r[Survey_Name_Season == all_survey_names[i],year]),max(distances_dissimilarities_allyears.r[Survey_Name_Season == all_survey_names[i],year])),
#       y = c(0.15,0.9)) +
  xlab("Year") +
ylab("beta-diversity") +
  facet_wrap(~Survey_Name_Season, ncol = 5) +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

print(points_wavylines_bray_total_year_reduced_gam_individual[[i]])

}

saveRDS(points_wavylines_bray_total_year_reduced_gam_individual, here::here("figures","points_wavylines_bray_total_year_reduced_gam_individual.Rds"))

#print to pdf
library(gridExtra)

ggsave(
   filename = here::here("figures","points_wavylines_bray_total_year_reduced_gam_individual.pdf"), 
   plot = marrangeGrob(points_wavylines_bray_total_year_reduced_gam_individual, nrow=1, ncol=1), 
   width = 8.5, height = 11
)

#Alternatively, split into 2 and use facet
#first 24
points_wavylines_bray_total_year_reduced_gam_individual_facet_1_24 <- ggplot() +
  geom_point(data = distances_dissimilarities_allyears.r[Survey_Name_Season  %in% all_survey_names[c(1:24)]],
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean), alpha = 0.7) +
    geom_line(data = year_survey_unit_expand.dt[Survey_Name_Season   %in% all_survey_names[c(1:24)]],
             aes(x = year,
                 y = bray_glm_mod_fit)) +
  geom_ribbon(data = year_survey_unit_expand.dt[Survey_Name_Season  %in% all_survey_names[c(1:24)]],
aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE), alpha=0.5) + #add standard error
    geom_smooth(data = distances_dissimilarities_allyears.r[Survey_Name_Season  %in% all_survey_names[c(1:24)]],
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 color = Survey_Name_Season,
                 fill = Survey_Name_Season), alpha = 0.4, method = "lm", linewidth = 0.5) +
    scale_color_manual(values =    linear_model_color_order_factorlevels$trend_color[1:24], name = "Survey Unit") +
    scale_fill_manual(values =    linear_model_color_order_factorlevels$trend_color[1:24], name = "Survey Unit") +
  theme_classic() +
  xlab("Year") +
ylab("β-diversity") +
  facet_wrap(~Survey_Name_Season, ncol = 4, scales = "free") +
  theme(axis.text = element_text(size = 8), axis.title = element_text(size = 12), legend.position = "null")

ggsave(points_wavylines_bray_total_year_reduced_gam_individual_facet_1_24, path =  here::here("figures"), filename = "points_wavylines_bray_total_year_reduced_gam_individual_facet_1_24.png", height = 11, width = 9)

#second half
distances_dissimilarities_allyears.r_second <- distances_dissimilarities_allyears.r[Survey_Name_Season  %in% all_survey_names[c(25:42)]]
distances_dissimilarities_allyears.r_second <- setorder(distances_dissimilarities_allyears.r_second,Survey_Name_Season)
year_survey_unit_expand.dt_second <- year_survey_unit_expand.dt[Survey_Name_Season   %in% all_survey_names[c(25:42)]]
year_survey_unit_expand.dt_second <- setorder(year_survey_unit_expand.dt_second, Survey_Name_Season)

points_wavylines_bray_total_year_reduced_gam_individual_facet_25_42 <- ggplot() +
  geom_point(data = distances_dissimilarities_allyears.r_second,
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean), alpha = 0.7) +
    geom_line(data = year_survey_unit_expand.dt_second,
             aes(x = year,
                 y = bray_glm_mod_fit)) +
  geom_ribbon(data = year_survey_unit_expand.dt_second,
aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE), alpha=0.5) + #add standard error
  geom_smooth(data = distances_dissimilarities_allyears.r_second,
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 color = Survey_Name_Season,
                 fill = Survey_Name_Season), alpha = 0.4, method = "lm", linewidth = 0.5) +
    scale_color_manual(values =    linear_model_color_order_factorlevels$trend_color[25:42], name = "Survey Unit") +
    scale_fill_manual(values =    linear_model_color_order_factorlevels$trend_color[25:42], name = "Survey Unit") +
  theme_classic() +
  xlab("Year") +
ylab("β-diversity") +
  facet_wrap(~Survey_Name_Season, ncol = 4, scales = "free") +
  theme(axis.text = element_text(size = 8), axis.title = element_text(size = 12))

ggsave(points_wavylines_bray_total_year_reduced_gam_individual_facet_25_42, path =  here::here("figures"), filename = "points_wavylines_bray_total_year_reduced_gam_individual_facet_25_42.png", height = 11, width = 9)

NA
NA

Region by region build up (for presentations)

#EBS ONLY (segue for chapter 3)
points_wavylines_bray_total_year_reduced_gam_colorbytrend_EBSonly <- ggplot() +
  geom_point(data = na.omit(distances_dissimilarities_allyears.r[survey_unit == "EBS"], cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), size = 2, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt[survey_unit == "EBS"], cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season)) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt[survey_unit == "EBS"], cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
    scale_color_manual(values =  "#cbbfde", name = "Survey Unit") +
  scale_fill_manual(values =  "#cbbfde", guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year])), y = c(0.49,0.95)) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam_colorbytrend_EBSonly

ggsave(points_wavylines_bray_total_year_reduced_gam_colorbytrend_EBSonly, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_colorbytrend_EBSonly.jpg", height = 6, width = 11, unit = "in")

#just one region that's homogenizing (SEUS-summer)
points_wavylines_bray_total_year_reduced_gam_colorbytrend_seus_summer <- ggplot() +
  geom_point(data = na.omit(distances_dissimilarities_allyears.r[survey_unit == "SEUS-summer"], cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), size = 2, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt[survey_unit == "SEUS-summer"], cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season)) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt[survey_unit == "SEUS-summer"], cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
    scale_color_manual(values =  "#e7ac5b", name = "Survey Unit") +
  scale_fill_manual(values =  "#e7ac5b", guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year])), y = c(0.49,0.95)) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam_colorbytrend_seus_summer

ggsave(points_wavylines_bray_total_year_reduced_gam_colorbytrend_seus_summer, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_colorbytrend_seus_summer.jpg", height = 6, width = 11, unit = "in")

#just one region that's differentiating and one that's homogenizing
points_wavylines_bray_total_year_reduced_gam_colorbytrend_NZ_WCSI_SEUS_summeronly <- ggplot() +
  geom_point(data = na.omit(distances_dissimilarities_allyears.r[survey_unit %in% c("NZ-WCSI","SEUS-summer")], cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), size = 2, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt[survey_unit %in% c("NZ-WCSI","SEUS-summer")], cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season)) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt[survey_unit %in% c("NZ-WCSI","SEUS-summer")], cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
    scale_color_manual(values =  c("#e7ac5b","#91c874"), name = "Survey Unit") +
  scale_fill_manual(values =  c("#e7ac5b","#91c874"), guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year])), y = c(0.49,0.95)) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam_colorbytrend_NZ_WCSI_SEUS_summeronly

ggsave(points_wavylines_bray_total_year_reduced_gam_colorbytrend_NZ_WCSI_SEUS_summeronly, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_colorbytrend_NZ_WCSI_SEUS_summeronly.jpg", height = 6, width = 11, unit = "in")

#just three regions, another one that's stable, (NEUS-Fall)
points_wavylines_bray_total_year_reduced_gam_colorbytrend_NZ_WCSI_SEUS_summer_NEUSfallonly <- ggplot() +
  geom_point(data = na.omit(distances_dissimilarities_allyears.r[survey_unit %in% c("NZ-WCSI","SEUS-summer", "NEUS-Fall")], cols = "year_adj"),
             aes(x = year,
                 y = bray_curtis_dissimilarity_total_mean,
                 fill = Survey_Name_Season), size = 2, shape = 21, color = "white") +
    geom_line(data = na.omit(year_survey_unit_expand.dt[survey_unit %in% c("NZ-WCSI","SEUS-summer", "NEUS-Fall")], cols = "year_adj"),
             aes(x = year,
                 y = bray_glm_mod_fit,
                 color = Survey_Name_Season)) +
  geom_ribbon(data = na.omit(year_survey_unit_expand.dt[survey_unit %in% c("NZ-WCSI","SEUS-summer", "NEUS-Fall")], cols = "year_adj"), aes(x = year, ymin=bray_glm_mod_fit-bray_glm_mod_fit_SE, ymax=bray_glm_mod_fit+bray_glm_mod_fit_SE, fill =  Survey_Name_Season), alpha=0.1) + #add standard error
    scale_color_manual(values =  c("#cbbfde","#e7ac5b", "#91c874"), name = "Survey Unit") +
  scale_fill_manual(values =  c("#cbbfde","#e7ac5b", "#91c874"), guide = "none") +
  theme_classic() +
  lims(x = c(min(distances_dissimilarities_allyears.r[,year]),max(distances_dissimilarities_allyears.r[,year])), y = c(0.49,0.95)) +
  xlab("Year") +
ylab("β-diversity") +
  theme(legend.position = "null", axis.text = element_text(size = 15), axis.title = element_text(size = 15))

points_wavylines_bray_total_year_reduced_gam_colorbytrend_NZ_WCSI_SEUS_summer_NEUSfallonly

ggsave(points_wavylines_bray_total_year_reduced_gam_colorbytrend_NZ_WCSI_SEUS_summer_NEUSfallonly, path = here::here("figures"), filename ="points_wavylines_bray_total_year_reduced_gam_colorbytrend_NZ_WCSI_SEUS_summer_NEUSfallonly.jpg", height = 6, width = 11, unit = "in")

NA
NA
NA

Merge BC versus Year plot with GAMS and Region vs. coefficient plot for LMs

START HERE

Is each survey better described by LM or by GAM?

# Function to compare GAM and linear model fits
compare_models <- function(subset_data) {
  # Fit GAM model with "fs" basis and single smooth term
  gam_model <- gam(bray_curtis_dissimilarity_total_mean ~ s(year, bs = "fs", m = 1), data = subset_data)
  
  # Fit linear model
  lm_model <- lm(bray_curtis_dissimilarity_total_mean ~ year, data = subset_data)
  
  # Calculate AICc values
  aicc_gam <- AICc(gam_model)
  aicc_lm <- AICc(lm_model)
  
  comparison <- ifelse((abs(aicc_gam-aicc_lm) < 2 | (aicc_gam == aicc_lm)),"Same",
                       ifelse(((abs(aicc_gam-aicc_lm) > 2 & (aicc_gam < aicc_lm))),"GAM",
                              ifelse(((abs(aicc_gam-aicc_lm) > 2 & (aicc_gam > aicc_lm))), "Linear",NA)))
  
   return(data.table(Survey_Name_Season = unique(subset_data$Survey_Name_Season)[1],
                    comparison = comparison,
                    aicc_gam = aicc_gam,
                    aicc_lm = aicc_lm))
}

GAM_LM_model_comparisons <- data.table()

# Compare models for each survey name and bind the results
for (i in 1:length(all_survey_names)) {
  GAM_LM_model_comparisons_single <- compare_models(distances_dissimilarities_allyears.r[Survey_Name_Season == all_survey_names[i]])
  GAM_LM_model_comparisons <- rbind(GAM_LM_model_comparisons, GAM_LM_model_comparisons_single)
}

Warning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NAWarning: Item 2 has 0 rows but longest item has 1; filled with NAWarning: Item 3 has 0 rows but longest item has 1; filled with NA

GAM_LM_model_comparisons

#round to 2 significant figures

#names of numeric columns
numeric_cols <- names(GAM_LM_model_comparisons)[sapply(GAM_LM_model_comparisons, is.numeric)]

# Apply signif() only to numeric columns
GAM_LM_model_comparisons[, (numeric_cols) := lapply(.SD, function(x) if (is.numeric(x)) signif(x, digits = 3) else x), .SDcols = numeric_cols]


View(GAM_LM_model_comparisons)

colnames(GAM_LM_model_comparisons) <- c("Survey","Model comparison", "GAM AICc","LM AICc")

fwrite(GAM_LM_model_comparisons, here::here("output","GAM_LM_model_comparisons.csv"))

SCRATCH OLD CODE These bar plots show LMER random slope and intercept and use alternate color schemes

BC_total_Dissimilarity_Coef_errorbar_reduced <- ggplot() +
    geom_errorbar(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, ymin = lwr, ymax = upr), fill = "grey", width = 0) + #add confidence intervals
  geom_point(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, size = years_sampled, fill = Directional_Change_sig, color = Directional_Change_sig), stat = 'identity', shape = 21) +
  scale_fill_manual(values = c("white","black","grey"), name = "Dissimilarity trend", guide="none") +
  scale_color_manual(values = c("black","black","grey"), name = "Dissimilarity trend", guide="none") +
  scale_size_binned(range = c(1,8), name = "Survey period length\n(years)") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = seq(-0.005, 0.0075, by = 0.0025), labels = c("-0.005","","0", "", "0.005",  "")) +
  xlab("Survey unit") +
  ylab("β-diversity trend") + #total BC dissimilarity trend
  coord_flip() +
  theme_classic() +
  theme(axis.text.y = element_text(colour = color_year_adj_order, face = "bold"), axis.title.y = element_blank(), axis.text.x = element_text(size = 15), axis.title.x = element_text(size = 15), legend.position = c(0.25,0.7), legend.direction = "vertical")

#pull legend for homogenization
directional_change_legend_plot <- BC_total_Dissimilarity_Coef_errorbar_reduced + 
  scale_fill_manual(values = c("white","black","grey"), name = "Dissimilarity trend") +
  scale_color_manual(values = c("black","black","grey"), name = "Dissimilarity trend") +
  scale_size_binned(range = c(1,8), name = "Survey period length", guide = "none") +
  theme(legend.position = "right", legend.background = element_rect(fill= "transparent")) +
  guides(colour = guide_legend(override.aes = list(size=6)))


BC_total_Dissimilarity_Coef_errorbar_reduced

#ALT grey scale
BC_total_Dissimilarity_Coef_errorbar_reduced_greyscale <- ggplot() +
    geom_errorbar(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, ymin = lwr, ymax = upr), fill = "grey", width = 0) + #add confidence intervals
  geom_point(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, size = years_sampled, fill = Directional_Change_sig, color = Directional_Change_sig), stat = 'identity', shape = 21) +
  scale_fill_manual(values = c("white","black","grey"), name = "Dissimilarity trend", guide="none") +
  scale_color_manual(values = c("black","black","grey"), name = "Dissimilarity trend", guide="none") +
  scale_size_binned(range = c(1,8), name = "Survey period length\n(years)") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = seq(-0.005, 0.0075, by = 0.0025), labels = c("-0.005","","0", "", "0.005",  "")) +
  xlab("Survey unit") +
  ylab("β-diversity trend") + #total BC dissimilarity trend
  coord_flip() +
  theme_classic() +
  theme(axis.text.y = element_text(face = "bold"), axis.title.y = element_blank(), axis.text.x = element_text(size = 15), axis.title.x = element_text(size = 15), legend.position = c(0.25,0.7), legend.direction = "vertical")

#"#73BA4D","#E0962C","#cbbfde"

BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend <- ggplot() +
    geom_errorbar(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, ymin = lwr, ymax = upr), fill = "grey", width = 0) + #add confidence intervals
  geom_point(data = model_coefs_reduced_length, aes(x = reorder(Survey_Name_Season, year_adj) , y = year_adj, label = Survey_Name_Season, size = years_sampled, color = Directional_Change_sig), stat = 'identity') +
  scale_color_manual(values = c("#73BA4D","#E0962C","#cbbfde"), name = "Dissimilarity trend", guide="none") +
  scale_size_binned(range = c(1,8), name = "Survey period length\n(years)") +
  geom_hline(yintercept = 0) +
  scale_y_continuous(breaks = seq(-0.005, 0.0075, by = 0.0025), labels = c("-0.005","","0", "", "0.005",  "")) +
  xlab("Survey unit") +
  ylab("β-diversity trend") + #total BC dissimilarity trend
  coord_flip() +
  theme_classic() +
  theme(axis.text.y = element_text(face = "bold"), axis.title.y = element_blank(), axis.text.x = element_text(size = 15), axis.title.x = element_text(size = 15), legend.position = c(0.3,0.8), legend.direction = "vertical", legend.text = element_text(size = 15), legend.title = element_text(size = 16))

directional_change_legend_plot_colorbytrend <- BC_total_Dissimilarity_Coef_errorbar_reduced_colorbytrend + 
  theme(legend.position = "right", legend.background = element_rect(fill= "transparent"), 
         legend.text = element_text(size = 15), legend.title = element_text(size = 16)) +
  guides(colour = guide_legend(override.aes = list(size=6)), size = "none")

Merge BC versus Year plot with GAMS and Region vs. coefficient plot for LMs


<!-- rnb-text-end -->


<!-- rnb-chunk-begin -->


<!-- rnb-source-begin eyJkYXRhIjoiYGBgclxuXG5CQ190b3RhbF9HQU1fTE1FUl9tZXJnZV9sZWdlbmQgPC0gZ2dkcmF3KHhsaW0gPSBjKDAsIDQwLjUpLCB5bGltID0gYygwLCAyMSkpICtcbiAgICBkcmF3X3Bsb3QocG9pbnRzX3dhdnlsaW5lc19icmF5X3RvdGFsX3llYXJfcmVkdWNlZF9nYW0sXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHggPSAxLCB5ID0gMSwgd2lkdGggPSAyMCwgaGVpZ2h0ID0gMjApICtcbiAgICBkcmF3X3Bsb3QoQkNfdG90YWxfRGlzc2ltaWxhcml0eV9Db2VmX2Vycm9yYmFyX3JlZHVjZWQgK1xuICAgICAgICB0aGVtZShsZWdlbmQua2V5LnNpemUgPSB1bml0KDAuNSwgJ2NtJyksICNjaGFuZ2UgbGVnZW5kIGtleSBzaXplXG4gICAgICAgIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfdGV4dChzaXplPTE2KSwgI2NoYW5nZSBsZWdlbmQgdGl0bGUgZm9udCBzaXplXG4gICAgICAgIGxlZ2VuZC50ZXh0ID0gZWxlbWVudF90ZXh0KHNpemU9MTQpKSwgI2NoYW5nZSBsZWdlbmQgdGV4dCBmb250IHNpemUpLFxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB4ID0gMjAsIHkgPSAxLCB3aWR0aCA9MTksIGhlaWdodCA9IDIwKSArXG4gICAgZHJhd19wbG90KGdldF9sZWdlbmQoZGlyZWN0aW9uYWxfY2hhbmdlX2xlZ2VuZF9wbG90ICsgXG4gICAgICB0aGVtZShsZWdlbmQua2V5LnNpemUgPSB1bml0KDAuNSwgJ2NtJyksICNjaGFuZ2UgbGVnZW5kIGtleSBzaXplXG4gICAgICAgIGxlZ2VuZC50aXRsZSA9IGVsZW1lbnRfdGV4dChzaXplPTE1KSwgI2NoYW5nZSBsZWdlbmQgdGl0bGUgZm9udCBzaXplXG4gICAgICAgIGxlZ2VuZC50ZXh0ID0gZWxlbWVudF90ZXh0KHNpemU9MTMpKSksICNjaGFuZ2UgbGVnZW5kIHRleHQgZm9udCBzaXplKVxuICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeCA9IDI3LCB5ID0gMTAsIHdpZHRoID0gMywgaGVpZ2h0ID0gMikgK1xuICBnZW9tX3RleHQoYWVzKHggPSAyLCB5ID0gMjAuNyksIGxhYmVsID0gKFwiYS5cIiksIHNpemUgPTgsIGZvbnRmYWNlID0gXCJib2xkXCIpICtcbiAgZ2VvbV90ZXh0KGFlcyh4ID0yMCwgeSA9IDIwLjcpLCBsYWJlbCA9IChcImIuXCIpLCBzaXplID04LCBmb250ZmFjZSA9IFwiYm9sZFwiKVxuXG5cbmdnc2F2ZShCQ190b3RhbF9HQU1fTE1FUl9tZXJnZV9sZWdlbmQsIHBhdGggPSBoZXJlOjpoZXJlKFwiZmlndXJlc1wiKSwgZmlsZW5hbWUgPSBcIkJDX3RvdGFsX0dBTV9MTUVSX21lcmdlX2xlZ2VuZC5wbmdcIiwgaGVpZ2h0ID0gOCwgd2lkdGggPSAxNCwgdW5pdHMgPSBcImluXCIpXG5cbiNBTFQgR1JFWSBTQ0FMRVxuQkNfdG90YWxfR0FNX0xNRVJfbWVyZ2VfbGVnZW5kX2dyZXlzY2FsZSA8LSBnZ2RyYXcoeGxpbSA9IGMoMCwgNDAuNSksIHlsaW0gPSBjKDAsIDIxKSkgK1xuICAgIGRyYXdfcGxvdChwb2ludHNfd2F2eWxpbmVzX2JyYXlfdG90YWxfeWVhcl9yZWR1Y2VkX2dhbV9ncmV5c2NhbGUsXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHggPSAxLCB5ID0gMSwgd2lkdGggPSAyMCwgaGVpZ2h0ID0gMjApICtcbiAgICBkcmF3X3Bsb3QoQkNfdG90YWxfRGlzc2ltaWxhcml0eV9Db2VmX2Vycm9yYmFyX3JlZHVjZWRfZ3JleXNjYWxlICtcbiAgICAgICAgdGhlbWUobGVnZW5kLmtleS5zaXplID0gdW5pdCgwLjUsICdjbScpLCAjY2hhbmdlIGxlZ2VuZCBrZXkgc2l6ZVxuICAgICAgICBsZWdlbmQudGl0bGUgPSBlbGVtZW50X3RleHQoc2l6ZT0xNiksICNjaGFuZ2UgbGVnZW5kIHRpdGxlIGZvbnQgc2l6ZVxuICAgICAgICBsZWdlbmQudGV4dCA9IGVsZW1lbnRfdGV4dChzaXplPTE0KSksICNjaGFuZ2UgbGVnZW5kIHRleHQgZm9udCBzaXplKSxcbiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgeCA9IDIwLCB5ID0gMSwgd2lkdGggPSAxOSwgaGVpZ2h0ID0gMjApICtcbiAgICBkcmF3X3Bsb3QoZ2V0X2xlZ2VuZChkaXJlY3Rpb25hbF9jaGFuZ2VfbGVnZW5kX3Bsb3QgKyBcbiAgICAgIHRoZW1lKGxlZ2VuZC5rZXkuc2l6ZSA9IHVuaXQoMC41LCAnY20nKSwgI2NoYW5nZSBsZWdlbmQga2V5IHNpemVcbiAgICAgICAgbGVnZW5kLnRpdGxlID0gZWxlbWVudF90ZXh0KHNpemU9MTUpLCAjY2hhbmdlIGxlZ2VuZCB0aXRsZSBmb250IHNpemVcbiAgICAgICAgbGVnZW5kLnRleHQgPSBlbGVtZW50X3RleHQoc2l6ZT0xMykpKSwgI2NoYW5nZSBsZWdlbmQgdGV4dCBmb250IHNpemUpXG4gICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHggPSAyNywgeSA9IDEwLCB3aWR0aCA9IDMsIGhlaWdodCA9IDIpICtcbiAgZ2VvbV90ZXh0KGFlcyh4ID0gMiwgeSA9IDIwLjcpLCBsYWJlbCA9IChcImEuXCIpLCBzaXplID04LCBmb250ZmFjZSA9IFwiYm9sZFwiKSArXG4gIGdlb21fdGV4dChhZXMoeCA9MjAsIHkgPSAyMC43KSwgbGFiZWwgPSAoXCJiLlwiKSwgc2l6ZSA9OCwgZm9udGZhY2UgPSBcImJvbGRcIilcblxuZ2dzYXZlKEJDX3RvdGFsX0dBTV9MTUVSX21lcmdlX2xlZ2VuZF9ncmV5c2NhbGUsIHBhdGggPSBoZXJlOjpoZXJlKFwiZmlndXJlc1wiKSwgZmlsZW5hbWUgPSBcIkJDX3RvdGFsX0dBTV9MTUVSX21lcmdlX2xlZ2VuZF9ncmV5c2NhbGUucG5nXCIsIGhlaWdodCA9IDgsIHdpZHRoID0gMTQsIHVuaXRzID0gXCJpblwiKVxuYGBgIn0= -->

```r

BC_total_GAM_LMER_merge_legend <- ggdraw(xlim = c(0, 40.5), ylim = c(0, 21)) +
    draw_plot(points_wavylines_bray_total_year_reduced_gam,
                                         x = 1, y = 1, width = 20, height = 20) +
    draw_plot(BC_total_Dissimilarity_Coef_errorbar_reduced +
        theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
        legend.title = element_text(size=16), #change legend title font size
        legend.text = element_text(size=14)), #change legend text font size),
                                         x = 20, y = 1, width =19, height = 20) +
    draw_plot(get_legend(directional_change_legend_plot + 
      theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
        legend.title = element_text(size=15), #change legend title font size
        legend.text = element_text(size=13))), #change legend text font size)
                                       x = 27, y = 10, width = 3, height = 2) +
  geom_text(aes(x = 2, y = 20.7), label = ("a."), size =8, fontface = "bold") +
  geom_text(aes(x =20, y = 20.7), label = ("b."), size =8, fontface = "bold")


ggsave(BC_total_GAM_LMER_merge_legend, path = here::here("figures"), filename = "BC_total_GAM_LMER_merge_legend.png", height = 8, width = 14, units = "in")

#ALT GREY SCALE
BC_total_GAM_LMER_merge_legend_greyscale <- ggdraw(xlim = c(0, 40.5), ylim = c(0, 21)) +
    draw_plot(points_wavylines_bray_total_year_reduced_gam_greyscale,
                                         x = 1, y = 1, width = 20, height = 20) +
    draw_plot(BC_total_Dissimilarity_Coef_errorbar_reduced_greyscale +
        theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
        legend.title = element_text(size=16), #change legend title font size
        legend.text = element_text(size=14)), #change legend text font size),
                                         x = 20, y = 1, width = 19, height = 20) +
    draw_plot(get_legend(directional_change_legend_plot + 
      theme(legend.key.size = unit(0.5, 'cm'), #change legend key size
        legend.title = element_text(size=15), #change legend title font size
        legend.text = element_text(size=13))), #change legend text font size)
                                x = 27, y = 10, width = 3, height = 2) +
  geom_text(aes(x = 2, y = 20.7), label = ("a."), size =8, fontface = "bold") +
  geom_text(aes(x =20, y = 20.7), label = ("b."), size =8, fontface = "bold")

ggsave(BC_total_GAM_LMER_merge_legend_greyscale, path = here::here("figures"), filename = "BC_total_GAM_LMER_merge_legend_greyscale.png", height = 8, width = 14, units = "in")